setup <- tempfile(fileext = ".r")
readLines(con = 'https://inkaverse.com/setup.r') |>
gsub("echo = FALSE", "echo = TRUE", x = _) |>
writeLines(setup)
source(setup)
library(gtsummary)
library(factoextra)
library(corrplot)Impact of water deficit on growth, productivity, and water use efficiency in potato genotypes (Solanum tuberosum L.)
1 Setup
2 Googlesheets connect
url <- "https://docs.google.com/spreadsheets/d/1dfgpmCKdPmxRHozrZp0iE_xMGsvKTIcztDpMWYSEGaY"
gs <- as_sheets_id(url)
# browseURL(url)3 Trait rename
geno <- c("CIP720088" = "G01"
, "CIP392797.22" = "G02"
, "CIP397077.16" = "G03"
, "CIP398192.213" = "G04"
, "CIP398180.612" = "G05"
, "CIP398208.704" = "G06"
, "CIP398098.119" = "G07"
, "CIP398190.89" = "G08"
, "CIP398192.592" = "G09"
, "CIP398201.510" = "G10"
, "CIP398203.244" = "G11"
, "CIP398203.5" = "G12"
, "CIP398208.219" = "G13"
, "CIP398208.33" = "G14"
, "CIP398208.620" = "G15")
traits <- c("SPAD_29" = "Chlorophyll concentration (SPAD) at 29 dap"
, "SPAD_59" = "Chlorophyll concentration (SPAD) at 59 dap"
, "SPAD_76" = "Chlorophyll concentration (SPAD) at 76 dap"
, "SPAD_83" = "Chlorophyll concentration (SPAD) at 83 dap"
, "HGT" = "Plant height (cm)"
, "RWC" = "Relative water content (%)"
, "LOP" = "Leaf osmotic potential (MPa)"
, "LDW" = "Leaf dry weight (g)"
, "SDW" = "Stem dry weight (g)"
, "RDW" = "Root dry weight (g)"
, "TDW" = "Tuber dry weight (g)"
, "NTUB" = "Tuber number (N°)"
, "TRS" = "Total transpiration (mL)"
, "LFA" = "Leaf area (cm^2^)"
, "RTL" = "Root length (cm)"
, "TDB" = "Total dry biomass (g)"
, "HI" = "Harvest index (HI)"
, "SLA" = "Specific leaf area (cm^2^/g)"
, "RCC" = "Relative chlorophyll content (RCC)"
, "WUE_B" = "Biomass water use efficiency (g/L)"
, "WUE_T" = "Tuber water use efficiency (g/L)"
)4 Import data
fb <- gs %>%
range_read("fb") %>%
rename_with("tolower") %>%
rename_with(~gsub("\\s+|\\.", "_", .)) %>%
mutate(treat = case_when(treat %in% "wellwater" ~ "WW"
, treat %in% "drydown" ~ "WD")) %>%
select(block, treat, genotype,
spad_29 = spad_29dap,
spad_59 = spad_59dap,
spad_76 = spad_76dap,
spad_83 = spad_83dap,
hgt = hgt_86dap,
rwc = rwc_84dap,
lop = op_84dap,
ldw = leafdw,
sdw = stemdw,
rdw = rootdw,
tdw = tubdw,
ntub,
trs = ttrns,
lfa = la,
rtl = rlg
) %>%
mutate(tdb = (ldw+sdw+rdw+tdw),
hi = tdw/(ldw+sdw+rdw+tdw),
sla = lfa/ldw,
rcc = (spad_83/lfa),
wue_b = (ldw+sdw+rdw+tdw)/trs,
wue_t = tdw/trs
) %>%
mutate(gnt = recode(genotype, !!!geno), .after = genotype) %>%
select(genotype, gnt, treat, block, everything()) %>%
mutate(across(where(is.character), as.factor))
# str(fb)
fb %>% web_table()5 Abbreviations
gs %>%
range_read("var") %>%
filter(include == "x") %>%
select(Variable, Abbreviation) %>%
kable()| Variable | Abbreviation |
|---|---|
| Soil Plant Analysis Development | SPAD |
| Plant height | HGT |
| Relative water content | RWC |
| Leaf osmotic potential | LOP |
| Leaf dry weight | LDW |
| Stem dry weight | SDW |
| Root dry weight | RDW |
| Tuber dry weight | TDW |
| Tuber number | NTUB |
| Root length | RTL |
| Total transpiration | TRS |
| Leaf area | LFA |
| Total dry biomass | TDB |
| Harvest index | HI |
| Specific leaf area | SLA |
| Water use efficiency | WUE |
| Tuber water use efficiency | TWUE |
6 Table 1
tab <- gs %>%
range_read("gnt") %>%
dplyr::select(Number, Genotypes, Adaptability, "Growing period" = GPL, "Heat tolerance" = Heat, "Dry matter (%)")
# tab %>% sheet_write(gs, data = ., sheet = "tab1")
tab %>% web_table()7 Table 2
tab <- fb %>%
select(!c(block, gnt, genotype)) %>%
tbl_summary(
by = treat,
statistic = list(all_continuous() ~ "{mean} ± {sd}"),
missing = "no") %>%
add_p() %>%
bold_labels() %>%
as_tibble() %>%
mutate(`**Characteristic**` = recode(
`**Characteristic**`
, "__spad_29__" = "Chlorophyll concentration (SPAD) at 29 dap"
, "__spad_59__" = "Chlorophyll concentration (SPAD) at 59 dap"
, "__spad_76__" = "Chlorophyll concentration (SPAD) at 76 dap"
, "__spad_83__" = "Chlorophyll concentration (SPAD) at 83 dap"
, "__hgt__" = "Plant height (cm)"
, "__rwc__" = "Relative water content (%)"
, "__lop__" = "Leaf osmotic potential (MPa)"
, "__ldw__" = "Leaf dry weight (g)"
, "__sdw__" = "Stem dry weight (g)"
, "__rdw__" = "Root dry weight (g)"
, "__tdw__" = "Tuber dry weight (g)"
, "__ntub__" = "Tuber number (N°)"
, "__trs__" = "Total transpiration (mL)"
, "__lfa__" = "Leaf area (cm^2^)"
, "__rtl__" = "Root length (cm)"
, "__tdb__" = "Total dry biomass (g)"
, "__hi__" = "Harvest index (HI)"
, "__sla__" = "Specific leaf area (cm^2^/g)"
, "__rcc__" = "Relative chlorophyll content (RCC)"
, "__wue_b__" = "Biomass water use efficiency (g/L)"
, "__wue_t__" = "Tuber water use efficiency (g/L)"
)) %>%
rename(Variable = `**Characteristic**`
, `Water deficit` = `**WD**, N = 75`
, `Well-Watered`= `**WW**, N = 75`
, `p-value` = `**p-value**`)
# tab %>% sheet_write(gs, data = ., sheet = "tab2")
tab %>% web_table()8 Figure 1
fts <- gs %>%
range_read("ftsw") %>%
filter(Treatment != "preharvest") %>%
tidyr::gather(key = day, value = fts, -ID, -Genotype, -Treatment)
model <- fts %>%
inti::mean_comparison(response = "fts"
, model_factors = "Treatment*day"
, c("Treatment", "day"))
plt1 <- model$comparison %>%
mutate(Treatment = case_when(
Treatment == "drydown" ~ "Water Deficit (WD)"
, Treatment == "wellwater" ~ "Well Watered (WW)"
)) %>%
mutate(across(day, as.numeric)) %>%
plot_smr(
type = "line", color = F
, x = "day"
, y = "fts"
, group = "Treatment"
, ylab = "Fraction of transpirable soil water"
, xlab = "Days"
, glab = "Treatment"
, legend = "none"
, ylimits = c(0, 1.08, 0.1)
, error = "ste"
) +
scale_x_continuous(breaks = c(0:100)*2
, limits = c(34, 88))
trns <- gs %>%
range_read("trans") %>%
filter(Treatment != "preharvest") %>%
tidyr::gather(key = day, value = trans, -ID, -Genotype, -Treatment) %>%
filter(day != "TOTAL") %>%
drop_na()
model <- trns %>%
inti::mean_comparison(response = "trans"
, model_factors = "Treatment*day"
, c("Treatment", "day"))
plt2 <- model$comparison %>%
mutate(Treatment = case_when(
Treatment == "drydown" ~ "Water Deficit (WD)"
, Treatment == "wellwater" ~ "Well Watered (WW)"
)) %>%
mutate(across(day, as.numeric)) %>%
plot_smr(type = "line", color = F
, x = "day"
, y = "trans"
, group = "Treatment"
, ylab = "Transpiration (mL)"
, xlab = "Days"
, glab = "Irrigation"
, ylimits = c(0, 700, 100)
, error = "ste"
) +
theme(axis.title.x=element_blank(),
axis.text.x=element_blank(),
axis.ticks.x=element_blank())
plot <- ggdraw(xlim = c(0, 0.5), ylim = c(0,0.9)) +
draw_plot(plt2, width = 0.5, height = 0.43 , x = 0.0, y = 0.47 ) +
draw_plot(plt1, width = 0.495, height = 0.5, x = 0.005, y = 0.0 ) +
draw_plot_label(
label = c("a", "b"),
x = c(0.465, 0.465),
y = c(0.79, 0.49))
ggsave2(plot = plot, "files/Fig1.tiff", width = 18, height = 15, units = "cm")
"files/Fig1.jpg" %>% include_graphics()9 Figure 2
# tdw
model <- fb %>%
inti::mean_comparison(response = "tdw"
, model_factors = "block + gnt*treat"
, c("gnt","treat"))
tdw <- model$comparison %>%
plot_smr(type = "bar"
, color = F
, x = "gnt"
, xlab = ""
, y = "tdw"
, ylab = "Tuber dry weight (g)"
, group = "treat"
, glab = "Treatment"
, legend = "none"
, error = "ste"
, ylimits = c(0, 100, 20)
)
# rcc
model <- fb %>%
inti::mean_comparison(response = "rcc"
, model_factors = "block + gnt*treat"
, c("gnt","treat"))
rcc <- model$comparison %>%
plot_smr(type = "bar"
, color = F
, x = "gnt"
, xlab = ""
, y = "rcc"
, ylab = "Relative chlorophyll content"
, ylimits = c(0, 0.1, 0.02)
, group = "treat"
, glab = "Treatment"
, legend = "none"
, error = "ste"
)
# hi
model <- fb %>%
inti::mean_comparison(response = "hi"
, model_factors = "block + gnt*treat"
, c("gnt","treat"))
hi <- model$comparison %>%
plot_smr(type = "bar"
, color = F
, x = "gnt"
, xlab = ""
, y = "hi"
, ylab = "Harvest Index"
, group = "treat"
, glab = "Treatment"
, legend = "none"
, ylimits = c(0, 1, 0.2)
, error = "ste"
)
# wue_t
model <- fb %>%
inti::mean_comparison(response = "wue_t"
, model_factors = "block + gnt*treat"
, c("gnt","treat"))
wue_t <- model$comparison %>%
plot_smr(type = "bar", color = F
, x = "gnt"
, xlab = ""
, y = "wue_t"
, ylab = "Tuber water use efficiency (g/L)"
, group = "treat"
, glab = "Treatment"
, legend = "none"
, ylimits = c(0, 10, 2)
, error = "ste"
)
# spad_29
model <- fb %>%
inti::mean_comparison(response = "spad_29"
, model_factors = "block + gnt*treat"
, c("gnt","treat"))
spad_29 <- model$comparison %>%
mutate(treat = case_when(
treat == "WD" ~ "Water Deficit (WD)"
, treat == "WW" ~ "Well Watered (WW)"
)) %>%
plot_smr(type = "line"
, x = "gnt"
, xlab = ""
, y = "spad_29"
, ylab = "SPAD at 29 dap"
, group = "treat"
, glab = "Treatment"
, legend = "top"
, ylimits = c(30, 70, 5)
, error = "ste"
, color = F
)
# spad_83
model <- fb %>%
inti::mean_comparison(response = "spad_83"
, model_factors = "block + gnt*treat"
, c("gnt","treat"))
spad_83 <- model$comparison %>%
plot_smr(type = "line", color = F
, x = "gnt"
, xlab = ""
, y = "spad_83"
, ylab = "SPAD at 83 dap"
, group = "treat"
, glab = "Treatment"
, legend= "none"
, ylimits = c(30, 70, 5)
, error = "ste"
)
plot <- ggdraw(xlim = c(0.0, 1), ylim = c(0.0, 1)) +
draw_plot(spad_29 , width = 0.49, height = 0.35, x = 0.0, y = 0.6) +
draw_plot(spad_83, width = 0.49, height = 0.3, x = 0.5, y = 0.6) +
draw_plot(rcc, width = 0.49, height = 0.3, x = 0.0, y = 0.3) +
draw_plot(tdw, width = 0.49, height = 0.3, x = 0.5, y = 0.3) +
draw_plot(hi , width = 0.49, height = 0.3, x = 0.0, y = 0.0) +
draw_plot(wue_t, width = 0.49, height = 0.3, x = 0.5, y = 0.0) +
draw_plot_label(
label = c("a", "b", "c", "d", "e", "f"),
x = c(0.06, 0.56, 0.06, 0.56, 0.06, 0.56),
y = c(0.89, 0.89, 0.59, 0.59, 0.29, 0.29))
ggsave(plot = plot, "files/Fig2.tiff", units = "cm", width = 28, height = 25)
"files/Fig2.jpg" %>% include_graphics()10 Figure 3
cr <- fb %>%
group_by(gnt, treat) %>%
summarise_all(list(mean), na.rm = TRUE) %>%
as.data.frame() %>%
select(!c( genotype, block)) %>%
unite(gnt, treat, col = "rn") %>%
column_to_rownames("rn") %>%
rename_with(toupper, where(is.numeric))
cp <- heatmaply::heatmaply_cor(
cor(cr),
k_col = 5,
dendrogram = c("column"),
grid_gap = 1,
cellnote = round(cor(cr),2),
cellnote_textposition = "middle center",
cellnote_size = 10,
file = "files/correlation.pdf"
)
fig <- magick::image_read_pdf("files/correlation.pdf") %>%
grid::rasterGrob() %>%
ggsave("files/Fig3.tiff", plot = ., width = 8.39, height = 5.26)
unlink("files/correlation.pdf")
"files/Fig3.jpg" %>% include_graphics()Where: Chlorophyll Concentration (SPAD), Plant height (HGT; cm), Relative water content (RWC; %), Leaf osmotic potential (LOP; MPa), Leaf dry weight (LDW; g), Stem dry weight (SDW; g), Root dry weight (RDW; g), Tuber dry weight (TDW; g), Tuber number (NTUB; N°), Total transpiration (TRS; mL), Leaf area (LFA; cm2), Root length (RTL; cm), Total dry biomass (TDB; g), Harvest index (HI), Specific leaf area (SLA; cm2g-1), Relative chlorophyll content (RCC), Biomass water use efficiency (WUEB; gL-1), Tuber water use efficiency (WUET; gL-1).
11 Figure 4
mvd <- fb %>%
select(!c(block, gnt)) %>%
group_by(treat, genotype) %>%
summarise(across(everything(), ~mean(.x, na.rm = T))) %>%
mutate(coln = paste(genotype, treat, sep = "_")) %>%
column_to_rownames("coln") %>%
select(!genotype) %>%
rename_with(toupper, where(is.numeric))
pca <- PCA(mvd, graph = F, scale.unit = TRUE, quali.sup = 1)
# PCA
tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=8*ppi, height=8*ppi, res=ppi)
plot.PCA(pca,
choix="var",
title="",
autoLab = "auto",
shadowtext = T)
graphics.off()
pt1 <- png::readPNG(tmp) %>%
grid::rasterGrob()
# sink("files/pca.txt")
# # Results
# summary(pca, nbelements = Inf, nb.dec = 2)
# # Correlation de dimensions
# dimdesc(pca)
# sink()
# Analysis de Hierarchical Clustering
clus <- HCPC(pca, nb.clust=-1, graph = F)
tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=8*ppi, height=8*ppi, res=ppi)
plot.HCPC(clus,
choice = "map",
legend = list(bty = "y", x = "topright"),
draw.tree = F)
graphics.off()
pt2 <- png::readPNG(tmp) %>%
grid::rasterGrob()
# sink("files/clu.txt")
# clus$call$t$tree
# clus$desc.ind
# clus$desc.var
# sink()
plot <- plot_grid(pt1, pt2, nrow = 1, labels = "auto")
ggsave2(plot = plot, "files/Fig4.tiff", units = "cm", width = 20, height = 11)
"files/Fig4.jpg" %>% include_graphics()Where: Chlorophyll Concentration (SPAD), Plant height (HGT; cm), Relative water content (RWC; %), Leaf osmotic potential (LOP; MPa), Leaf dry weight (LDW; g), Stem dry weight (SDW; g), Root dry weight (RDW; g), Tuber dry weight (TDW; g), Tuber number (NTUB; N°), Total transpiration (TRS; mL), Leaf area (LFA; cm2), Root length (RTL; cm), Total dry biomass (TDB; g), Harvest index (HI), Specific leaf area (SLA; cm2g-1), Relative chlorophyll content (RCC), Biomass water use efficiency (WUEB; gL-1), Tuber water use efficiency (WUET; gL-1).
12 Supplementary Table 1
vars <- names(fb)
model <- 5:length(vars) %>% map(function(x) {
fb %>%
H2cal(trait = vars[x]
, gen.name = "genotype"
, rep.n = 4
, fixed.model = "0 + (1|block) + genotype"
, random.model = "1 + (1|block) + (1|genotype)"
, emmeans = T
, plot_diag = T
, outliers.rm = T
)
})
tabsmr <- 1:length(model) %>% map(function(x) {
model[[x]][["tabsmr"]]
}) %>%
bind_rows() %>%
mutate(trait = recode(trait, !!!traits)) %>%
select(trait, mean, std, min, max, V.g, V.e, repeatability) %>%
mutate(across(where(is.numeric), ~ round(.x, digits = 2))) %>%
mutate(trait = toupper(trait)) %>%
mutate(trait = recode(trait, !!!traits)) %>%
rename(Trait = trait)
tabsmr %>% web_table()
# tabsmr %>% sheet_write(gs, data = ., sheet = "tabS1")13 Supplementary Figure 1
var <- get_pca_var(pca)
tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=3.8*ppi, height=10*ppi, res=ppi)
corrplot(var$cor,
method="number",
tl.col="black",
tl.srt=45,)
graphics.off()
pt1 <- png::readPNG(tmp) %>%
grid::rasterGrob(interpolate = TRUE)
pt2 <- fviz_eig(pca,
addlabels=TRUE,
hjust = 0.05,
barfill="white",
barcolor ="darkblue",
linecolor ="red") +
ylim(0, 50) +
labs(
title = "PCA - percentage of explained variances",
y = "Variance (%)") +
theme_minimal()
pt3 <- fviz_contrib(pca,
choice = "var",
axes = 1,
top = 10,
fill="white",
color ="darkblue",
sort.val = "desc") +
ylim(0, 18) +
labs(title = "Dim 1 - variables contribution")
pt4 <- fviz_contrib(pca,
choice = "var",
axes = 2,
top = 10,
fill="white",
color ="darkblue",
sort.val = "desc") +
ylim(0, 18) +
labs(title = "Dim 2 - variables contribution")
plot <- ggdraw(xlim = c(0.0, 1.0), ylim = c(0, 1.0))+
draw_plot(pt1, width = 0.4, height = 0.99, x = 0.6, y = 0.0) +
draw_plot(pt2, width = 0.6, height = 0.34, x = 0.03, y = 0.66) +
draw_plot(pt3, width = 0.6, height = 0.34, x = 0.03, y = 0.33) +
draw_plot(pt4, width = 0.6, height = 0.34, x = 0.03, y = 0.0) +
draw_plot_label(
label = c("a", "b", "c", "d"),
x = c(0.005, 0.005, 0.005, 0.65),
y = c(0.999, 0.67, 0.34, 0.999))
ggsave2(plot = plot, "files/SFig1.tiff", height = 25, width = 30, units = "cm")
"files/SFig1.jpg" %>% include_graphics()14 Supplementary Figure 2
fig <- magick::image_read("files/potatos.png") %>%
grid::rasterGrob() %>%
ggsave("files/SFig2.tiff", plot = ., width = 15, height = 15)
"files/SFig2.jpg" %>% include_graphics()15 Convert figures to jpg format
library(tidyverse)
figs <- "files/" %>% list.files(pattern = "tiff", full.names = T)
convert <- 1:length(figs) %>% map(\(x) {
figname <- figs[x] %>% basename() %>% gsub("tiff", "jpg", .)
figs[x] %>%
tiff::readTIFF(., native=TRUE) %>%
jpeg::writeJPEG(., target = file.path("files", figname))
})